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Understanding the amazingly complex human cerebral cortex requires a map (or parcellation) of its major subdivisions, 
known as cortical areas. Making an accurate areal map has been a century-old objective in neuroscience. Using multi- 
modal magnetic resonance images from the Human Connectome Project (HCP) and an objective semi-automated 
neuroanatomical approach, we delineated 180 areas per hemisphere bounded by sharp changes in cortical architecture, 
function, connectivity, and/or topography in a precisely aligned group average of 210 healthy young adults. We 
characterized 97 new areas and 83 areas previously reported using post-mortem microscopy or other specialized study- 
specific approaches. To enable automated delineation and identification of these areas in new HCP subjects and in future 
studies, we trained a machine-learning classifier to recognize the multi-modal ‘fingerprint’ of each cortical area. This 
classifier detected the presence of 96.6% of the cortical areas in new subjects, replicated the group parcellation, and 
could correctly locate areas in individuals with atypical parcellations. The freely available parcellation and classifier 
will enable substantially improved neuroanatomical precision for studies of the structural and functional organization 
of human cerebral cortex and its variation across individuals and in development, aging, and disease. 


Neuroscientists have long sought to subdivide the human brain into a 
mosaic of anatomically and functionally distinct, spatially contiguous 
areas (cortical areas and subcortical nuclei), as a prerequisite 
for understanding how the brain works. Areas differ from their 
neighbours in microstructural architecture, functional specialization, 
connectivity with other areas, and/or orderly intra-area topographic 
organization (for example, the map of visual space in visual cortical 
areas)!~, Accurate parcellation provides a map of where we are in 
the brain, enabling efficient comparison of results across studies and 
communication among investigators; as a foundation for illuminating 
the functional and structural organization of the brain; and as a means 
to reduce data complexity while improving statistical sensitivity and 
power for many neuroimaging studies. 

The human cerebral cortex has been estimated to contain any- 
where from ~50 (ref. 1) to ~200 (refs 3, 4) areas per hemisphere. 
However, attaining a consensus whole-cortex parcellation has been 
difficult because of practical and technical challenges that we address 
here. 

Most previous parcellations were based on only one neurobiological 
property (such as architecture, function, connectivity or topography), 
and many cover only part of the cortex. Combining multiple properties 
provides complementary as well as confirmatory information, as 
different properties distinguish different sets of areal boundaries, 
and more confidence can be placed in boundaries that are consistent 
across multiple independent properties. We analysed all four proper- 
ties across all of neocortex in both hemispheres, using new or refined 
methods applied to the uniquely rich repository of exceptionally 
high-quality magnetic resonance imaging (MRI) data provided by 
the Human Connectome Project (HCP), which benefited from major 
advances in image acquisition and preprocessing ®. Architectural 
measures of relative cortical myelin content and cortical thickness were 


derived from T1-weighted (T1w) and T2-weighted (T2w) structural 
images*?"°. Cortical function was measured using task functional MRI 
(tfMRI) contrasts from seven tasks!". Resting-state functional MRI 
(rfMRI) revealed functional connectivity of entire cortical areas plus 
topographic organization within some areas. 

Previous parcellations typically used either fully automated 
algorithmic approaches, or else manual or partly automated 
neuroanatomical approaches in which neuroanatomists delineated 
areal borders, documented areal properties, and identified areas after 
consulting prior literature. Here we combined both approaches. For 
the initial parcellation, we adapted a successful observer-independent 
semi-automated neuroanatomical approach for generating post- 
mortem architectonic parcellations'*!? to non-invasive neuroimaging. 
We used an algorithm to delineate potential areal borders (transitions 
in two or more of the cortical properties described above), which two 
neuroanatomists (authors M.EG. and D.C.V.E.) then interpreted, 
documenting areal properties and identifying areas relative to the 
extant neuroanatomical literature. We then used a fully automated 
algorithmic approach, training a machine-learning classifier to 
delineate and identify cortical areas in individual subjects based 
on multi-modal areal fingerprints, allowing the parcellation to be 
replicated in new subjects and studies. 

Prior parcellations have either used small numbers of individuals 
or group averages that are ‘blurry’ from inaccurate alignment of brain 
areas across subjects. We aligned cortical data using ‘areal features, 
including maps of relative myelin content and resting state networks 
that are more closely tied to cortical areas than are the folding patterns 
typically used for alignment'*. The markedly improved intersubject 
cortical alignment using cortical folding, myelin, and resting state fMRI 
enabled us to generate the ‘typical subject’s’ parcellation from a highly 
detailed 210-subject group average data set. 
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Group map reproducibility of fine details 

We analysed two independent groups of HCP subjects—210P 
(‘parcellatiom) and 210V (‘validation’)—aligned using areal-feature- 
based registration (called “MSMAII, see Methods section on image 
preprocessing). Figure 1 illustrates the consistency of fine spatial 
patterns for maps reflecting relative myelin content (left panels) and 
a task-f{MRI language-related activation (right panels). The maps 
are strikingly similar across the 210P and 210V groups, including 
variations in relative myelin content within the primary somatosen- 
sory cortex related to somatotopic organization (white and black arrows 
in left panels, see legend) and small features in the task fMRI maps 
(white ellipses in right panels). Supplementary Results and Discussion 
1.1 and associated Supplementary Figs 1-5 include more examples 
of such cross-group consistency for architecture (myelin and cortical 
thickness), function (tfMRI contrast maps), and two resting state 
connectivity measures. 

Because of the areal-feature-based alignment, maps of average 
cortical folding in this study are much blurrier than are maps of areal 
properties because folding patterns and area locations are imperfectly 
correlated*'*"4 (for example, compare Supplementary Fig. 7e and 7j 
(group average and individual subject folding) in Supplementary 
Results and Discussion 1.3). Group average folding patterns remain 
sharp mainly in early sensory areas, where areal locations and folds 
are tightly correlated (for example, central and calcarine sulci, 
see Supplementary Fig. 1, rows 3 and 4, in Supplementary Results and 
Discussion 1.1). The regional difference in sharpness between maps of 
areal properties and folding highlights the importance of alignment 
based on areal features, rather than folding patterns, as a prerequisite for 
accurately parcellating group average data. The high spatial resolution 
of the HCP’s MRI images and lack of aggressive spatial smoothing? 
prior to group averaging also contribute to making our maps 
substantially sharper than those from traditional neuroimaging studies. 

Quantitatively, the 210P and 210V group average datasets were highly 
correlated across the cortical surface (r= 0.998 for myelin; r=0.994 for 
cortical thickness after correction for folding-related effects; r= 0.996 
and r=0.979 for two folding-related measures (FreeSurfer’s ‘sule and 
‘curv, respectively); r=0.995, r=0.984 and r=0.944 for the maximum, 
median and minimum, respectively, of the task fMRI contrasts, and a 
median reproducibility of r= 0.989 for two measures of resting state 
connectivity). These excellent map reproducibilities provide confidence 
that the parcellation will reflect the areal pattern of typical subjects in 
the healthy young adult population. See Methods section on modalities 
for parcellation, and Supplementary Methods 1.3-3.4 for the methods 
used to generate these maps. 


Myelin map 


Figure 1 | Consistency of fine spatial details in independent group 
averages. Relative myelin content maps (left hemisphere) and task 
fMRI contrast beta maps from the LANGUAGE story contrast (right 
hemisphere) on inflated (columns 1 and 3) and flattened surfaces 
(columns 2 and 4). Rows 1 and 2 are the group averages of the 210P and 
210V data sets, respectively. White and black arrows indicate consistent 
variations in myelin content within primary somatosensory cortex that 
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A 180-area group average parcellation 

To identify transitions representing candidate areal boundaries, we 
designed and implemented a semi-automated, quantitative approach 
adapted for multi-modal neuroimaging data represented on two- 
dimensional cortical surface models (see Methods section on the 
gradient-based parcellation approach and Supplementary Methods 
4.1-5.3). The approach is similar in spirit to a highly successful 
semi-automated observer-independent approach!*"!>. However, 
instead of objectively identifying potential areal borders in postmortem 
histological sections, we identified them algorithmically on the cortical 
surface by computing the first derivative of each areal feature map (its 
spatial gradient magnitude)'®. Candidate borders were then interpreted 
by the neuroanatomists to exclude artefacts. Each area's properties were 
documented (in the Supplementary Neuroanatomical Results), and 
putative areas were related to the extant neuroanatomical literature. 

These semi-automated approaches contrast with classical 
observer-dependent parcellation approaches!” that have relied on 
visual inspection to locate often subtle transitions in cortical architec- 
ture and with some modern observer-dependent retinotopic parcella- 
tion methods!”'®, They also differ from fully automated, unsupervised 
methods’? *! in which the outcomes depend heavily on algorithmic 
input parameters (for example, thresholds or number of requested 
clusters) and are not validated by a neuroanatomist. 

Area 55b illustrates our multi-modal gradient-based parcellation 
approach using gradients of three areal feature maps (see Fig. 2). Area 
55b is a small, elongated, and notably distinct area (outlined in black 
or white) bounded by the frontal eye field (FEF) and premotor eye field 
(PEF), primary motor cortex (4), ventral premotor cortex (6v), and pre- 
frontal areas 8Av and 8C. In the myelin map (Fig. 2a), area 55b is lightly 
myelinated and lies between moderately myelinated areas FEF (above) 
and PEF (below), just anterior to heavily myelinated primary motor 
cortex (area 4). Thus, area 55b is surrounded on three sides by myelin 
gradients (Fig. 2e). Area 55b is strongly activated in the ‘Story versus 
Baseline’ task contrast from the HCP’s ‘LANGUAGE task (Fig. 2b) 
and is entirely surrounded by a strong gradient for this task contrast 
(Fig. 2f). It also has distinctive functional connectivity, as revealed by a 
seed location (lightly myelinated area PSL) selectively connected with 
55b (Fig. 2c) and a different seed location (heavily myelinated area 
LIPv) strongly connected with FEF and PEF (Fig. 2d) but not with 55b. 
The result is strong mean gradients in dense functional connectivity 
surrounding 55b (Fig. 2g). Ref. 22 illustrated area 55b on a schematic 
surface map (Fig. 2h) as a lightly myelinated area bounded on three 
sides by more heavily myelinated areas. Because of the similarity to the 
dorsal portion of 55b in ref. 22, we use the same name. 


Task fMRI story 
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are correlated with somatotopy (see Supplementary Neuroanatomical 
Results 6 and Supplementary Neuroanatomical Results Fig. 8). The 

white oval indicates a small, sharp, and reproducible feature in the right 
hemisphere of the LANGUAGE story contrast. Relative myelin content 
will hereafter be referred to as myelin (see legend of Supplementary Fig. 1 
in Supplementary Results and Discussion 1.1). Data at http://balsa.wustl. 
edu/WDpX. 
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Functional connectivity maps 
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Figure 2 | Parcellation of exemplar area 55b using multi-modal 
information. The border of 55b is indicated by a white or black outline. 
a, Myelin map. b, Group average beta map from the LANGUAGE Story 
versus Baseline task contrast. c, d, Functional connectivity correlation 
maps from a seed in area PSL (white sphere, arrow) (c) and a seed in 
area LIPv (white sphere, arrow) (d). e, Gradient magnitude of the myelin 
map shown in a. f, Gradient magnitude of the LANGUAGE Story versus 


To generate the complete parcellation of 180 areas and area 
complexes in each hemisphere, we adopted a systematic, objective, and 
quantitative approach (see the gradient-based parcellation approach 
section in the Methods and in Supplementary Methods 5.1-5.3). Our 
major criteria, met in nearly all cases, included: (i) spatially overlapping 
gradient ‘ridges’ between each pair of areas for at least two independent 
areal feature maps; (ii) similar gradient ridges present in roughly cor- 
responding locations in both hemispheres; (iii) gradients that were not 
correlated with artefacts; and (iv) robust and statistically significant 
cross-border differences in the feature maps. Another consideration 
(but not a requirement) was whether published evidence exists for a 
boundary in an approximately corresponding location. Studies with 
publicly available parcellations registered onto atlas surfaces* were 
directly compared with our data; however, most regions required 
indirect comparisons with published figures (for example, Fig. 1h). 


Baseline task contrast shown in b. g, Mean gradient magnitude of the 
functional connectivity dense connectome (see section on modalities for 
parcellation in the Methods). h, A dorsal schematic view of the prefrontal 
cortex as parcellated in ref. 22, in which shading indicates the amount 

of myelin found using histological stains of cortical grey matter. Data at 
http://balsa.wustl.edu/Qv4P. 


Initial areal boundaries meeting these criteria were delineated by two 
neuroanatomists (authors M.EG. and D.C.V.E.). 

In a second computational stage, the path of each manually drawn 
border was optimized algorithmically using gradients of the most 
informative feature maps selected by the neuroanatomists (those with 
visually obvious gradients and differences across the border). These 
feature maps were confirmed to have robust and statistically significant 
differences across the final border. The semi-automated gradient-based 
parcellation approach is further described in Supplementary Methods 
5.1-5.3), and the entire semi-automated process is illustrated for area 
V1 in Supplementary Neuroanatomical Results 1; other sections of this 
document describe and illustrate the information used to delineate and 
the literature used to name all 180 cortical areas. 

Figure 3 shows the multi-modal cortical parcellation in the left 
and right hemispheres on inflated and flattened surfaces, with areal 


The HCP’s multi-modal cortical parcellation (HCP_MMP1.0) 
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Figure 3 | The HCP’s multi-modal parcellation, version 1.0 (HCP_ 
MMP1.0). The 180 areas delineated and identified in both left and right 
hemispheres are displayed on inflated and flattened cortical surfaces. Black 
outlines indicate areal borders. Colours indicate the extent to which the 
areas are associated in the resting state with auditory (red), somatosensory 


positive 
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(green), visual (blue), task positive (towards white), or task negative 
(towards black) groups of areas (see Supplementary Methods 5.4). 
The legend on the bottom right illustrates the 3D colour space used 
in the figure. Data at http://balsa.wustl.edu/WN56. 
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boundaries delineated by black contours. A total of 180 areas and areal 
complexes per hemisphere is near the higher end of earlier estimates 
noted above**. We consider 180 likely to be a lower bound, as some 
parcels are probably complexes of multiple areas (for example, based 
on finer-grained published parcellations, and other regions that 
suffer from reduced sensitivity due to fMRI signal loss). Many areas 
(83/180) were assigned names based on published parcellations from 
dozens of separate studies that used a variety of invasive or specialized 
methods (see Table 2 in Supplementary Neuroanatomical Results), 
reflecting how far the field has been from a consensus neuroanatomical 
parcellation. Some of the newly described 97 areas have de novo names 
(for example, DVT for the dorsal visual transitional area), while others 
represent finer-grained parcellations of previously reported areas 
(for example, area 31 into areas 31a, 31pd, and 31pv). A few repre- 
sent complexes in which a published finer grained parcellation was 
not visible in our data (for example, areas 29 and 30 combined into 
area the retrosplenial complex (RSC)), but these may be again sub- 
divided once higher resolution data is available. The 180 areas differ 
widely in their shapes, sizes, and the positions of their borders relative 
to cortical folds. 

The parcellation in Fig. 3 is coloured to reflect each area’s degree of 
association in the resting state (determined using multiple regression, 
see Supplementary Methods 5.4) with five functionally specialized 
groups of areas: early auditory (red), early somatosensory/motor 
(green), and early visual areas (blue). These represent the three 
dominant input streams to the brain. Also used were two core groups 
of cognitive areas that are strongly anti-correlated in our data, the task 
positive network (towards white) and task negative (also called the 
default mode) network (towards black). Hence, the strongly bluish, 
greenish, and reddish regions are predominantly but not exclusively 
associated with visual, somatosensory-motor, and auditory processing, 
respectively. Qualitatively, the predominantly unimodal regions appear 
to collectively occupy less than half of the neocortical sheet. Areas 
probably more strongly biomodal include blue-green areas such as LIPv 
and MT (visual and somatosensory-motor) and purple areas such as 
POS2 and RSC (visual and auditory). The remaining regions form a 
complex mosaic, with some intermixing of lighter (task-positive) and 
darker (task-negative) areas along with many lighter or darker pastel 
hues suggestive of ‘cognitive’ areas that may be preferentially associated 
with one or another sensory modality. The bilateral symmetry of func- 
tional organization is striking, in that nearly all areas have qualitatively 
very similar hues in the left and right hemispheres. However, interesting 
colour asymmetries occur in a few areas, especially language-related 
areas 55b, PSL, SFL, and 44 and their right hemisphere homologues, 
which also have asymmetric task-fMRI functional profiles 
(see Supplementary Neuroanatomical Results 8, 15, 21 and 22). 

Internal heterogeneity is evident in some cortical areas, particularly 
those with topographically organized representations. In the 
somatosensory-motor strip (largely architecturally defined soma- 
tosensory and motor areas 3a, 3b, 1, 2, and 4), we identified five clearly 
defined topographic subareas in resting state and task {MRI data 
(see Supplementary Neuroanatomical Results 6 and the associated 
Supplementary Fig. 8). In this parcellation we treat topographic 
subdivisions as ‘subareas’ rather than calling them full ‘areas. For 
visual cortex, its visuotopic organization revealed a set of hemifield 
representations in each hemisphere, something not achieved in 
previous unsupervised resting state functional connectivity-based 
parcellations!?7!°, Also, ultrahigh-field MRI reveals sub-areal cortical 
organization along both laminar**”°, and columnar?®?’ axes, so our 
parcellation represents one of many important levels of granularity in 
brain organization. 


Cross -validation of the parcellation 

The initial statistical analysis used in the semi-automated parcella- 
tion was circular, to the extent that the 210P dataset was used for both 
creating and testing the parcellation. Hence, we carried out an additional 
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statistical cross-validation using the 210V dataset and a comprehensive 
set of feature maps (see the statistical cross validation of the multi-modal 
parcellation section of the Methods and Supplementary Results and 
Discussion 1.2). This analysis also reveals which areal properties were 
most useful in defining areal boundaries (a condensed representa- 
tion of the detailed information provided in the Supplementary 
Neuroanatomical Results). Supplementary Fig. 6 in Supplementary 
Results and Discussion 1.2 shows four independent categories of 
features: cortical thickness, myelin maps, task fMRI, and resting state 
fMRI, and how many of these categories showed robust and statistically 
significant differences across each areal border. Fully 96% of areal 
borders had robust effect sizes (Cohen’s d > 1) in two or more feature 
categories and all were statistically significant after correcting for 
multiple comparisons in two or more feature categories in cross-border, 
across-subject t-tests. Resting state fMRI was the most useful category, 
followed by task fMRI, myelin maps, and lastly cortical thickness, 
which was consistent with the neuroanatomists’ observations and 
documentation in the Supplementary Neuroanatomical Results. 


Exemplar parcellation-based analyses 

Spatial smoothing is often used to increase the signal-to-noise ratio 
(SNR) in neuroimaging analyses, to try to compensate for inaccurate 
registration of brain areas, and/or to satisfy statistical assumptions. 
However, smoothing blurs data across boundaries between areas 
(on the surface) and tissue compartments (in the volume). An areal 
parcellation enables area-wise analyses (averaging data within each 
area), thereby improving SNR and statistical power without the 
deleterious effects of spatial smoothing (to the extent that properties 
within an area are uniform). Parcellation dramatically reduces data 
dimensionality, illustrated here using the HCP’s myelin, thickness, task, 
and resting state data (Fig. 4). 

The ‘dense’ (vertex-wise) myelin map shown in Fig. 4a has ~30,000 
surface grey matter vertices per hemisphere, whereas a ‘parcellated’ 
myelin map (Fig. 4f) shows the same overall pattern with 180 cortical 
areas (vertices within an area have the same value, see also Fig. 4g 
for parcellated cortical thickness). Example dense and parcellated 
task fMRI analysis contrast maps (Figs. 4b, c, LANGUAGE Story 
versus Baseline) can be represented as a single column (white) in a 
180-area by 86-task-contrast matrix (Fig. 4d). Parcellated analyses 
hold great promise for task fMRI studies, as they improve the signal- 
to-noise ratio by averaging fMRI time series within parcels prior to 
fitting the task design, increasing Z statistics (Fig. 4e). Parcellation is 
effectively a neurobiologically constrained smoothing approach that 
also increases statistical power by efficiently consolidating otherwise 
non-independent statistical tests. This approach will benefit studies 
aimed at understanding the functional and structural organization 
of the brain in health or disease at an area-wise level (studies that 
currently summarize results using three-dimensional coordinates in 
a standardized stereotaxic space). Parcellated analyses also aid in the 
clarity and efficiency of communicating results (for example: “area 55b 
in the left hemisphere showed a statistically significant +1% BOLD 
activation in my language task”). 

Parcellated analyses are comparably useful when characterizing 
structural or functional connectivity, as previously recognized*>®. 
Preprocessing of HCP data results in fMRI data represented as 
‘grayordinates’ (cortical grey matter surface vertices and subcortical 
grey matter voxels”). A dense connectome, containing connectivity 
between all pairs of 91,282 grayordinates is ~3.3 x 10°-fold larger 
than an area-wise parcellated connectome for ~500 areas (connectivity 
between all pairs of areas), yet the parcellated connectome captures 
the neurobiologically relevant variance at the areal level. Parcellated 
connectomes are illustrated using a seed location in area PGi (black 
dot) for full correlation (Fig. 4h) and partial correlation (Fig. 4i) 
functional connectivity brain maps together with their associated 
parcellated connectome matrices (Fig 4j, full correlation below 
and partial correlation above the diagonal). In both cases, the task 


© 2016 Macmillan Publishers Limited, part of Springer Nature. All rights reserved. 


Dense myelin Dense task 


WAX X 


A / x 
20 0 90 
= 


g Parcellated thickness h 


4% 
Figure 4 | Example parcellated analyses using the HCP’s multi-modal 
cortical parcellation. a, Dense myelin maps on lateral (top) and medial 
(bottom) views of inflated left hemisphere. b, c, Example dense (b) and 
parcellated (c) task fMRI analysis (LANGUAGE story versus baseline) 
expressed as Z statistic values. d, The entire HCP task fMRI battery's 
Z statistics for 86 contrasts (47 unique, see section on modalities for 
parcellation in the Methods) analysed in parcellated form and displayed 
as a matrix (rows are parcels, columns are contrasts, white outline 
indicates the map in c). e, A major improvement in Z statistics from fitting 
task designs on parcellated time series instead of fitting them on dense 
time series and then parcellating afterwards (blue points are 


negative (default mode) network is evident, though the partial 
correlation connectome is much sparser than the full correlation 
connectome. 


Individuals with atypical areal patterns 

The precisely aligned group average multi-modal cortical parcellation 
represents the overall spatial arrangement of cortical areas in the 
‘typical’ individual from a healthy young adult population. However, 
we found atypical topological arrangements of some areas in some 
individuals that are discernible across multiple modalities, including 
resting-state networks, task-fMRI activations, and myelin maps. 
Distinguishing genuinely atypical areal topologies from inadequately 
aligned typical patterns depended on the MSMAIlareal-feature-based 
registration to align cortical areas precisely. We summarize key findings 
here and extensively characterize this important phenomenon in the 
Supplementary Results and Discussion 1.3. 

Previously described area 55b and neighbouring areas FEF and 
PEF showed particularly notable individual differences in topological 
arrangements. For the 210P subjects, 89% showed the typical 
configuration in the left hemisphere (area 55b bordered by area PEF 
inferiorly and area FEF superiorly, as in Fig. 2), which was well aligned 
with group average area 55b after MSMAII registration. However, in one 
subgroup (4%, n =9), a patch having the multi-modal characteristics 
of area 55b is shifted superiorly relative to the upper limb subregion 
of sensori-motor cortex (Supplementary Figure 7 in Supplementary 
Results and Discussion 1.3). In another subgroup (6%, n= 12), area 
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Parcellated task Parcellated vs dense task 
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360 parcels x 86 task contrasts; note the upward tilting deviation from the 
red line). f, Parcellated myelin maps. g, A parcellated folding-corrected 
cortical thickness map (in mm). h, i, Parcellated functional connectivity 
maps on the brain (seeded from area PGi, black dot). These parcellated 
connectomes are computed using either full or partial correlation (see 
Supplementary Methods 7.1). In both cases, the task negative (default 
mode) network is apparent. j, A parcellated connectome matrix view 
with the full correlation connectome below and the partial correlation 
connectome above the diagonal (white line shows the displayed partial 
correlation brain map). Data at http://balsa.wustl.edu/RGOx. 


55b is split into two pieces by a merger of areas FEF and PEE rather 
than the typical splitting of FEF and PEF by 55b (Supplementary 
Fig. 8 in Supplementary Results and Discussion 1.3). Such topological 
deviations in individual subjects’ areal maps raise intriguing questions 
for future exploration. They also cannot be corrected by a topology- 
preserving registration aimed at aligning individual subjects’ areas with 
the group average ‘atlas’ parcellation. Thus, we introduce an alternative 
fully automated cortical parcellation approach that can identify and 
delineate both typical and atypical areas in individual subjects that were 
nota part of the original 210P group. 


Automated individual-subject parcellation 

The semi-automated neuroanatomical approach described above is 
impractical for de novo individual subject parcellation of all ~1,100 
HCP subjects having complete MRI datasets so as to identify the 
atypical areal topologies mentioned above. Instead, we developed an 
automated method for generating individual subject parcellations 
based on a supervised machine learning classifier previously used to 
identify resting state functional networks in individual subjects”. In 
our case, the areal classifier learns the multi-modal ‘areal fingerprint 
of each cortical area that distinguishes it from surrounding cor- 
tex. Based on multi-modal feature maps that represent the areal 
properties of architecture, function, connectivity, and topography, 
the areal classifier returns a prediction (0% to 100%) that each area 
exists at a given cortical surface vertex. The highest prediction value 
across areas at each vertex is used to generate the individual subject 
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parcellation (see the cortical areal classifier section in Methods and 
in Supplementary Methods 6.1-6.8). Once trained using the 210P 
subjects (and a separate ‘29T” group of test subjects, see the subjects and 
acquisitions section of the Methods), the areal classifier should be 
able to use only the multi-modal areal fingerprints that it has learned 
to reproduce the parcellation in an independent group of validation 
subjects (210V). 

A critical early test of the areal classifier was whether it could 
accurately and reliably map areas that are not aligned with the 
population-based atlas parcellation after MSMAII areal-feature-based 
alignment (see Supplementary Results and Discussion 1.4). Examples 
of successful classification of areas 55b, FEF, and PEF are shown in 
Supplementary Fig. 9 of the Supplementary Results and Discussion 
for typical subjects, shifted 55b subjects, and split 55b subjects. In each 
illustrated case, the classifier correctly identified 55b and its neighbours 
(as assessed by the neuroanatomists’ inspection of the multi-modal 
areal features shown in the figure). Supplementary Fig. 10 in the 
Supplementary Results and Discussion 1.4 shows that these atypical 
55b topologies and classifications are stable across widely spaced repeat 
scanning sessions in a ‘test-retest’ group of 27 subjects (see Methods 
section on subjects and acquisition). 


Areal detection and parcellation consistency 

Another critical test of both the parcellation and the areal classifier 
is the classifier’s performance in detecting the 180 cortical areas in 
individual subjects, particularly in independent validation subjects 
that were not used to generate the parcellation or train the classifier. 
The top two rows of Fig. 5 show the performance of the classifier in 
detecting each area (see the cortical areal classifier section in Methods). 
Importantly, the classifier aims to detect whole areas based on their 
multi-modal fingerprints, rather than detecting differences in areal 
features across paired areal boundaries as was done in the cross- 
validation analysis (Supplementary Fig. 6 in Supplementary Results 
and Discussion 1.2). The overall areal detection rate was 98.0% of 
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Figure 5 | Areal detection rates, probabilistic areas, and parcellation 
reproducibility. Rows 1 (210P) and 2 (210V) show the individual subject 
areal detection rates (see Methods section on cortical areal classifier) 

as parcellated maps. Most areas are yellow (100%), and the minimum 
detection rate across both rows was 73%. Rows 3 and 4 illustrate 
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all areas across all subjects for the 210P parcellation and training 
dataset (row 1) and 96.6% for the independent 210V validation 
dataset (row 2), indicating excellent overall performance of the areal 
classifier. 

The areal classifier was used to generate probabilistic maps of 
each cortical area (illustrating residual variability in spatial location 
after MSMAII areal feature-based registration), and to assess the 
reproducibility of the parcellation in the independent 210V dataset. 
Rows 3 and 4 of Fig. 5 show strikingly similar probabilistic maps of 
8 non-overlapping areas with differing degrees of spatial variability 
(V1, 4, RSC, MT, LIPv, TEla, 46, and 10r) from the 210P and 210V 
groups. All probability maps were combined to produce a group 
maximum probability map (MPM), where the area with the highest 
probability at each vertex was found. Row 5 shows the original semi- 
automated parcellation borders, and row 6 compares the group MPM 
maps from 210P (blue) and 210V (red), with purple representing 
overlapping vertices. The borders in Row 6 are almost entirely purple, 
indicating very high reproducibility of the group MPM maps (r=0.965, 
Dice = 0.960, see the cortical areal classifier section of Methods). This 
reproducibility is similar to that of the original group average feature 
maps discussed above. The correlation of the original semi-automated 
parcellation (row 5) with the 210P group MPM (row 6) was r= 0.913, 
Dice = 0.902, indicating that the classifier made modest adjustments 
to better fit the data. We predict there will be very high reproducibility 
of the parcellation across the rest of the ~1,100 subject HCP dataset. 
Example individual subject parcellations and their reproducibility 
based on repeated scan sessions are shown in Supplementary Fig. 11 of 
Supplementary Results and Discussion 1.4. The individual parcellations 
are reasonably reproducible (median r= 0.77, Dice = 0.72) but, 
unsurprisingly, not as reproducible as the group parcellations, which 
benefit from averaging across many subjects. Other analyses yield 
interesting information about the sizes of cortical areas in the group 
average and variability in areal size across individuals (Supplementary 
Results and Discussion 1.5). 


probabilistic maps of areas V1, 4, RSC, MT, LIPv, TEla, 46, and 10r for 
the 210P (row 3) and 210V (row 4) groups. Row 5 shows the original 
parcellation derived from the semi-automated neuroanatomical approach. 
Row 6 shows the group MPM maps from 210P (blue), 210V (red), and 
their overlap (purple). Data at http://balsa.wustl.edu/WL8m. 
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Generalizing the classifier for future studies 

In contrast to the semi-automated approach (above) where neuroanat- 
omists chose the information to delineate and identify the 180 cortical 
areas in group-average data (see Supplementary Neuroanatomical 
Results), the areal classifier automatically determines (without human 
intervention) what information is most useful for delineating and 
identifying these cortical areas in individual subjects. As illustrated in 
Supplementary Figure 12 of Supplementary Results and Discussion 1.6, 
the areal classifier uses the task fMRI data least, perhaps because task 
fMRI feature maps are noisier in individual subjects than other feature 
maps and their information content is largely redundant with the 
resting state data”. This finding is important for generalizability of 
the areal classifier to other studies because replicating the custom- 
ized, hour-long HCP task fMRI battery is unlikely to be feasible for 
most neuroimaging studies. Ideally, the areal classifier would be able 
to perform nearly as well relying only on architecture, connectivity, 
and topography. Accordingly, we trained the classifier again on the 
210P dataset, but omitted the task {MRI-based feature maps. When 
trained this way, the classifier indeed performed nearly as well as 
when all features were used, detecting 97.6% of areas in 210P (versus 
98.0% using all features) and 96.4% of areas in 210V (versus 96.6% 
using all features). Hence, we anticipate that the areal classifier will 
generalize to other studies that acquire the following core set of MRI 
images: high-resolution Tlw and T2w; spin echo-based b0 field map; 
and extensive fMRI data acquired using ‘multiband’ pulse sequences to 
improve spatial and temporal resolution’ (see Supplementary Results 
and Discussion 2.3 and 2.9). These are the same image acquisition 
requirements as the HCP’s minimal preprocessing pipelines’ and the 
MSMaAll areal feature-based registration pipeline'* (Supplementary 
Methods 2.4). Future studies adhering to these image acquisition guide- 
lines will be able to use the unified framework of the HCP’s analysis 
pipelines to automatically generate individualized parcellated analy- 
ses from unprocessed MRI images, a major advance over traditional 
neuroimaging methods that have often relied on comparisons with 
Brodmann’ hand drawn parcellation published in 1909 (ref. 1). 


Discussion 
We have produced a population-based 180-area per hemisphere human 
cortical parcellation using exceptionally high quality multimodal data 
from hundreds of Human Connectome Project subjects aligned using 
an improved areal feature-based cross-subject alignment method 
(MSMAI)). Inspired by an observer-independent post-mortem architec- 
tural parcellation approach!’, we developed a semi-automated neuro- 
anatomical approach adapted to non-invasively acquired multi-modal 
MRI data. Although algorithms determined the final areal borders, 
the multi-modal data were carefully interpreted by neuroanatomists, 
the properties of each cortical area were documented, and each area 
was named in relation to the extant neuroanatomical literature (see 
Supplementary Neuroanatomical Results). A cross-validation showed 
that the areas forming the parcellation were robustly and statistically 
significantly different from their neighbours across multiple modalities. 
We identify this parcellation as HCP-MMP1.0 (Human Connectome 
Project Multi-Modal Parcellation version 1.0), making the version 1.0 
designation because we anticipate future refinements as better data 
become available (see Supplementary Results and Discussion 2.1). 
Unexpectedly, we discovered that despite improved intersubject 
alignment, some areas have atypical topological arrangements in 
some subjects, which we demonstrated for areas 55b, FEF, and PEF. 
We developed a fully automated method for parcellating individual 
subjects based on a machine learning classifier that can cope with 
this kind of individual variability. The areal classifier detected 96.6% 
of individual subject cortical areas in new subjects, including atypical 
areas, and replicated the group parcellation in an independent sample. 
Though we made extensive use of the HCP’s specialized task fMRI 
battery when generating the parcellation, we showed that task {MRI 
data is not essential for future studies aiming to use the areal classifier 
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to automatically define the cortical areas in their subjects. Instead, it 
suffices to acquire the same core set of MRI images needed for the rest 
of the HCP’s software pipelines. 

By generating a robust neuroanatomical map of human neocortical 
areas—a century-old aim of neuroscience—and providing methods 
for mapping these areas in any individual undergoing study with non- 
invasive neuroimaging, the present work represents a major advance 
relative to previous human cortical parcellations. The overall approach 
described here shows that we can produce sharp, reproducible brain 
images across multiple non-invasive neuroimaging modalities. We can 
generate a highly reproducible and generalizable cortical parcellation 
through state-of-the-art methods of data acquisition, preprocessing, 
and analysis designed to compensate for individual variability and 
thereby minimize blurring of images. These improvements, together 
with the new parcellation, make it desirable to use spatial localiza- 
tion methods that move beyond the traditional use of stereotaxic 
coordinates combined with Brodmann areal assignments to charac- 
terize centers of cortical activation in fMRI studies. From a neuro- 
anatomical perspective, there has often been substantial uncertainty 
whether any two neuroimaging studies have found results in the same 
cortical areas or not. The situation is analogous to astronomy in which 
ground-based telescopes produced relatively blurry images of the sky 
before the advent of adaptive optics and space telescopes. 

Many topics are discussed further in the Supplementary Results 
and Discussion 2.1-2.10 (for example, avenues for improving the 
parcellation and other issues left for future work, further discussion of 
the neuroscientific implications of our results, and additional datasets 
that could profitably be linked to our parcellation). As the topographic 
organization of higher cognitive areas becomes better understood, some 
parcels currently considered to be full areas may later be considered 
to be subareas of larger topographically organized cortical areas 
(analogous to somatotopic subregions of topographically organized 
sensory and motor areas illustrated in Supplementary Neuroanatomical 
Results 6). Though our use of multiple modalities probably mitigates 
this issue relative to traditional uni-modal parcellations, the extent to 
which the human multi-modal cortical parcellation may be revised 
along such lines remains a question for future work using the state 
of the art methods mentioned above (see Supplementary Results and 
Discussion 2.8). 

The MSMAII registration and the areal classifier are or will soon 
be freely available on GitHub; the visualization tool Connectome 
Workbench is on http://humanconnectome.org; and the parcellation, 
data, and scenes for reproducing each of the figures are in the BALSA 
database?!. These tools provide a neuroanatomical foundation, enabling 
the identification of cortical areas when reporting results or thinking 
about and discussing brain organization in relation to studies of human 
cognition, lifespan, and disease. Several additional interesting avenues 
of investigation are now open. The ability to discriminate individual 
differences in the location, size, and topology of cortical areas from dif- 
ferences in their activity or connectivity should facilitate the dissection 
of how each property is related to behaviour and genetic underpinnings, 
for example, in learning disabilities or those with distinctive cognitive 
traits. The ability to non-invasively and automatically delineate cortical 
areas in living subjects may have clinical implications, for example by 
providing neurosurgeons with detailed, individualized maps of the 
brains on which they operate. There are also important implications for 
our understanding of human cortical evolution. The dramatic expansion 
in neocortex along the human lineage occurred mainly in higher cogni- 
tive regions of lateral prefrontal, parietal, and temporal cortices™!3"?"?. 
Comparisons with nonhuman primates, including marmosets and 
macaques (both widely used in invasive studies), and great apes, may 
yield new insights regarding the emergence of new cortical areas and 
the divergences in areal functions, which collectively led to the cognitive 
capabilities that make us uniquely human as a species and as individuals. 

Note added in proof: A related paper on the neuroimaging approach 
used by the Human Connectome Project may be found in ref. 34. 
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In addition, we note that FreeSurfer uses an algorithm to label gyri and 
sulci automatically in individual subjects based on manually generated 
training labels?" that is similar in spirit to our areal classifier. Also, 
the FreeSurfer surface modelling noted in the Methods draws from 
methods summarized in ref. 36. 


Online Content Methods, along with any additional Extended Data display items and 
Source Data, are available in the online version of the paper; references unique to 
these sections appear only in the online paper. 
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METHODS 


Subjects and acquisition. A total of 449 young adult twins and non-twin siblings 
(ages 22-35) from the Human Connectome Project (HCP) were scanned according 
to the HCP’s acquisition protocol’. The MRI acquisition included collecting 
Tlw and T2w structural images, task-based and resting state-based fMRI images, 
diffusion-weighted images, and b0 field maps. Images were acquired at high spatial 
and temporal resolution on a customized Siemens 3 tesla (3T) scanner and with 
customized slice accelerated sequences for fMRI (see Supplementary Methods 
1.1-1.2). All subjects from the HCP 500-subject data release (July, 2014) having 
complete fMRI sessions were included. They were divided into two independent 
groups of 210 subjects that shared no family members between them, together 
with a remaining group of 29 test (29T) subjects that shared family members with 
210P but not 210V. The first group of subjects (210P, 130 females, 80 males) was 
used for creating the parcellation and training the areal classifier, which also made 
use of the 29T group to avoid overfitting. The second group of subjects (210V, 116 
females, 94 males) was used only for statistical cross-validation of the parcellation, 
areal classifier detection rates in independent subjects, and group parcellation 
reproducibility measures. A test-retest group of 27 subjects scanned twice through 
the entire MRI protocol and independently processed through the HCP pipelines 
was used for individual subject reproducibility measures. Subject recruitment 
procedures and informed consent forms, including consent to share de-identified 
data, were approved by the Washington University institutional review board. 
Datasets were de-identified and are publicly shared on the ConnectomeDB 
database (https://db.humanconnectome.org). 

Image preprocessing. Spatial image preprocessing (distortion correction and 
image alignment) was carried out using the HCP’s spatial minimal preprocessing 
pipelines’. These pipelines maximize alignment across image modalities, mini- 
mize distortions relative to the subject’s anatomical space, and minimize spatial 
smoothing (blurring) of the data. The data were projected into the 2mm standard 
CIFTI grayordinates space, which includes cortical grey matter surface vertices 
and subcortical grey matter voxels’. This offers substantial improvements in spa- 
tial localization over traditional volume-based analyses, enabling more accurate 
cross-subject and cross-study registrations and avoiding smoothing that mixes 
signals across differing tissue types or between nearby cortical folds. Additionally, 
we did minimal smoothing within the CIFTI grayordinates space to avoid mixing 
across areal borders prior to parcellation. 

For cross-subject registration of the cerebral cortex, we used a two-stage process 
based on the multimodal surface matching (MSM) algorithm! (see Supplementary 
Methods 2.1-2.5). An initial ‘gentle’ stage, constrained only by cortical folding 
patterns (FreeSurfer’s ‘sule measure), was used to obtain approximate geographic 
alignment without overfitting the registration to folding patterns, which are not 
strongly correlated with cortical areas in many regions. Previously, we found that 
more aggressive folding-based registration (either MSM-based or FreeSurfer- 
based) slightly decreased cross-subject task-fMRI statistics, suggesting that 
aligning cortical folds too tightly actually reduces alignment of cortical areas!*. A 
second, more aggressive stage used cortical areal features to bring areas into better 
alignment across subjects while avoiding neurobiologically implausible distortions 
or overfitting to noise in the data. The areal features used were myelin maps, resting 
state network maps computed with weighted regression (an improvement over 
dual regression?” described in the Supplementary Methods 2.3) and resting state 
visuotopic maps (see Supplementary Methods 4.4). Areal distortion was measured 
by taking the log base-2 of the ratio of the registered spherical surface tile areas to 
the original spherical surface tile areas. The mean (across space) of the absolute 
value of the areal distortion averaged across subjects from both registration 
stages was 30% less than the standard FreeSurfer folding-based registration and 
the maximum (across space) of this measure was 54% less. Despite less overall 
distortion, the areal-feature-based registration delivers substantially more accurate 
registration of cortical areas than does FreeSurfer folding-based registration as 
judged by cross-subject task {MRI statistics, an areal feature that was not used 
to drive the registration'*. Because MSM registration preserves topology and is 
relatively gentle (it does not tear or distort the cortical surface in neurobiologi- 
cally implausible ways), it is unable to align some cortical areas in some subjects 
where the areal arrangement differs from the group average (see Supplementary 
Results and Discussion 1.3-1.4 for more details on atypical areas). Group average 
registration drift away from the gentle folding-based geographic alignment was 
removed from the surface registration’? (see Supplementary Methods 2.5) to enable 
comparisons of this dataset with datasets registered using different areal features 
(for example, retinotopically defined areas). Group average registration drift is any 
consistent effect of the registration during template generation on the mean size, 
shape, or position of areas on the sphere (as opposed to the desired reductions in 
cross-subject variation). An obvious example is the 37% increase in average brain 
volume produced by registration to MNI space’. Uncorrected drifts during surface 
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template generation can cause apparent changes in cortical areal size, shape, and 
position when comparing across studies. 

Resting state fMRI data were denoised for spatially specific temporal artefacts 
(for example, subject movement, cardiac pulsation, and scanner artefacts) using 
the ICA+FIX approach, which includes detrending the data and aggressively 
regressing out 24 movement parameters*””. We avoided regressing out the ‘global 
signal’ (mean grey-matter time course) from our data because preliminary analyses 
showed that this step shifted putative connectivity-based areal boundaries so that 
they lined up less well with other modalities, likely because of the strong areal 
specificity of the residual global signal after ICA+FIX clean up. Task fMRI data 
were temporally filtered using a high pass filter. More details on resting state and 
task fMRI temporal preprocessing are described in the Supplementary Methods 
1.6-1.8. Substantial spatial smoothing was avoided for both datasets, and all images 
were intensity normalized to account for the receive coil sensitivity field. Artefact 
maps of large vein effects, fMRI gradient echo signal loss, and surface curvature 
were computed as described in Supplementary Methods 1.9. 

Modalities for parcellation. The multi-modal cortical parcellation used informa- 
tion related to the four areal properties of architecture, function, connectivity, and 
topography”. Architecture was measured using T1w/T2w myelin content maps plus 
cortical thickness maps with surface curvature regressed out*® (Supplementary 
Methods 1.5). Function was measured using task-fMRI responses to 7 tasks in 
86 task contrasts (47 unique; 39 were sign-reversed contrasts). Effect size maps 
(beta maps) after correction for the receive field were used instead of Z statistic 
maps because we were interested in regional differences in the magnitude of 
the BOLD (blood oxygen level dependent) signal change induced by the tasks, 
rather than differences in the significance of the BOLD signal change. Functional 
connectivity was measured using pairwise Pearson correlation of the denoised 
resting state time series of each pair of grayordinates. Topographic organization 
was explored using resting state time series in visual cortex, with spatial regressors 
representing polar angle and eccentricity patterns in area V1 combined with a 
modified ‘dual-regression-like’ approach that weights each surface vertex according 
to the cortical surface area that it represents (see Supplementary Methods 4.4). The 
semi-automated multi-modal parcellation was generated using group average data 
for all of these modalities from the 210P group of subjects (see Supplementary 
Methods 3.1-3.3 for details on how the group averages were created for each 
modality). The reproducibility of these group average maps was assessed by 
correlating the spatial maps for the 210P and 210V groups (see Supplementary 
Results and Discussion 1.1). 

The gradient-based parcellation approach. Classically, cortical areas have 
been defined based on sharp changes in one or more of the areal properties of 
architecture, function, connectivity, and topography. Traditionally, this relied 
heavily on visual inspection, until more objective and quantitative approaches 
became available!>!3, One highly successful approach to post-mortem architectural 
parcellation involves computing a dissimilarity metric, (the Mahalanobis distance) 
between neighbouring feature profiles generated from segmented histological 
images and testing for statistically significant and large spikes in dissimilarity that 
indicate putative areal boundaries. For in vivo data, a similarly powerful approach 
involves taking the first derivative (the spatial gradient) of a measure of interest 
along cortical surface and using the gradient magnitude to objectively identify 
locations where the measure is changing rapidly. One can then draw putative areal 
boundaries along the resulting gradient ridges!”'®'°. Here we combined elements 
of both approaches in a multi-modal context to generate semi-automatically 
drawn areal borders that were then evaluated statistically. Gradients were 
computed for architectural, functional, connectivity, and topographic modalities 
(see Supplementary Methods 4.1-4.4). 

To incorporate expert knowledge and priors from the neuroanatomical 
literature into the parcellation process, the neuroanatomists (authors M.E.G and 
D.C.V.E.) evaluated the multi-modal neuroimaging data and its gradients to define 
initial areal borders based on the following criteria. (1) Presence of a co-localized 
gradient ridge in at least two independent modalities was taken as strong evidence 
of an areal border, and the vast majority of areal borders satisfied this criterion. 
(2) Presence of corresponding gradients in the left and right hemispheres provided 
further evidence for a genuine areal border. For the vast majority of borders, 
the same modalities yielded robust gradients in both hemispheres. We did not 
find strong evidence for an area present in one hemisphere that was absent in 
the other (though a few areas show hemispheric asymmetries in their functional 
‘signature’ and/or in their spatial relationships with neighbouring areas). (3) We 
ignored gradients clearly attributable to imaging artefacts (see Supplementary 
Neuroanatomical Results for details). (4) Cortex on opposite sides of the border 
needed to differ robustly and significantly in the areal features used to delineate the 
border. (5) Confidence was increased if prior literature described a corresponding 
areal border. (6) Early runs of a supervised machine learning algorithm (see the 
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cortical areal classifier section of the Methods below) needed to be able to learn to 
distinguish each cortical area from its neighbours in a large majority of individual 
subjects based on individual subject multi-modal features (the early runs were only 
done using 210P and 29T, keeping 210V independent for later analyses). After the 
neuroanatomists delineated the initial areal borders and chose the important areal 
features that defined them, an automated algorithm then optimized the border 
placement so that it followed the most probable path based on the chosen areal 
features (see Supplementary Methods 5.1-5.3). The Supplementary 
Neuroanatomical Results documents the information that was used to distinguish 
each of the 180 areas from its neighbours. 

The neuroanatomists named areas based on previous parcellations whenever 

a reasonable match to the literature could be made. In some cases, areal identi- 
fication was based on the similarity of the area's properties relative to previously 
reported areas (for example, area 4, primary motor cortex, is known to be heavily 
myelinated and thick; area V2 has a mirror-image visuotopic map relative to 
neighbouring area V1). In most cases, however, the information used to describe 
previous cortical areas (for example, cytoarchitecture) was not available in the HCP 
data, and areal identification mainly reflected spatial correspondences relative to 
cortical folding patterns (if reliable for that region of cortex) or spatial relationships 
between neighbouring cortical areas. The strongest evidence for areal identifica- 
tion came from studies that provided surface-based probabilistic or maximum 
probability maps, ideally also registered using areal features and dedrifting of 
templates**. In these cases, we directly compared these data with our data and 
show the degree of overlap in the Supplementary Neuroanatomical Results. When 
such data were unavailable, we used published information to the degree feasible 
(see Supplementary Methods 5.3 for limitations of non-surface-based/not publicly 
available data) to make areal identifications or to describe new areas that had not 
previously been identified. The information used to name each cortical area is 
described in the Supplementary Neuroanatomical Results. 
Statistical cross validation of the multi-modal parcellation. Once the parcellation 
has been created, parcellated representations of data from each modality can be 
generated using either the group parcellation or the individual subject parcellations. 
For the statistical cross-validation, we created parcellated myelin, cortical thickness, 
task fMRI, and resting state functional connectivity datasets using the semi- 
automated multimodal group parcellation (see Supplementary Methods 7.1). For 
myelin and cortical thickness, we simply averaged the values of the dense individual 
subject maps within each area. For task fMRI, we averaged the time series within 
each area prior to computing task statistics (to benefit from the SNR improvements 
of parcellation demonstrated in Fig. 4e). For the same reason, we averaged resting 
state time series within each parcel prior to computing functional connectivity to 
form a parcellated functional connectome. 

For each pair of areas that shared a border in the parcellation, we computed a 
paired samples two-tailed t-test across subjects on these parcellated data for each 
feature (ignoring tests that involved the diagonal in the resting state parcellated 
functional connectome). We thresholded these tests at the Bonferroni-corrected 
significance level of P< 9 x 10 * (number of area pairs across both hemispheres 
(1,050) x number of features (266) x number of tails (2) x 0.05) and an effect size 
threshold of Cohen's d > 1. We grouped the features into 4 independent catego- 
ries (cortical thickness, myelin, task fMRI, and resting state fMRI) to determine 
for each area pair whether it showed robust and statistically significant 
differences across multiple modalities. For more details, see Supplementary 
Methods 7.2. 

The cortical areal classifier. We used a supervised machine learning classifier to 
automatically delineate and identify each cortical area from its neighbours across 
a large majority of individual subjects based on multi-modal information. Besides 
validating the robustness of the parcellation, this provides useful information about 
each individual subject’s parcellation, along with an approach to generalizing the 
parcellation to other datasets. To automatically parcellate individual subjects, we 
adapted the multi-layer perceptron used by ref. 29 to delineate and identify seven 
resting state networks more accurately than simpler linear methods including 
dual regression. We used the multi-layer perceptron to classify all 180 areas in our 
parcellation using multi-modal feature maps and relied on two neuroanatomically 
sensible assumptions to simplify the problem. (1) After areal feature-based registra- 
tion (MSMAII), we assumed that each cortical area was approximately in the same 
general location across subjects (for example, we don’t expect to find V1 outside 
the occipital lobe). This also means that we consider widely separated regions 
having similar multi-modal areal fingerprints to be distinct cortical areas even if 
they have similar architecture, coactivation in functional tasks, and belong to the 
same resting state network. These assumptions allowed us to reduce the overall 
classification problem to a set of 180 classification problems per hemisphere, each 
involving discrimination of one area from the areas around it. (2) Also, instead 
of classifying each area from all of its neighbours specifically (one class for the 


area plus one class for each neighbouring area), we set up the problem as a binary 
classification (the most robust kind of classification problem), classifying each 
area from all of the surrounding cortex as a single alternate class. This surrounding 
cortex represents a ‘searchlight’ for the area, and this searchlight was the group 
parcel location plus a 30 mm radius surrounding the group parcel in all directions 
across the surface (meaning that for a 10 mm circular area, the searchlight would be 
a circle of 70mm in diameter, still a quite large region of cortex). The 30 mm radius 
(geodesic distance computed on the group average mid-thickness surface corrected 
for vertex area loss due to averaging) was chosen because it easily encompassed 
the individual variation in area 55b in the 210P group (55b approaches a worst 
case because it is a relatively small and highly variable cortical area). The training 
labels were the group area from the semi-automated parcellation (class 1), and the 
remaining cortex in the searchlight (class 2). 

The features used by the classifier covered the same set of modalities used for 
the original parcellation, including architectural measures of myelin and cortical 
thickness with curvature regressed out; task fMRI maps (redundant information 
was reduced and SNR increased with a d= 20 ICA-decomposition run on the 
task contrast beta maps, see Supplementary Methods 6.4); the 77 surface-related 
resting state fMRI network maps computed on individual subjects using weighted 
regression from an overall d= 137 group ICA; five visuotopic topographic maps 
transformed into a format interpretable by the classifier; and maps of artefacts 
that the classifier used to interpret differences in areal features due to artefac- 
tual effects (see Supplementary Methods 6.3-6.5 for further description of each 
modality’s classifier features). These 112 multi-modal feature maps were gener- 
ated for each vertex in each of the 449 subjects and the 27 repeated subjects, with 
each hemisphere processed separately. Other than the 30 mm radius searchlight 
region of interest (ROI), the classifier has no spatial concept of where the area 
should be (it operates independently on each vertex and only knows what the 
area's fingerprint looks like in the feature space). Consequently, special consid- 
eration was given to the spatial visuotopic patterns, which were transformed into 
maps whose values reflected the alternating mirror symmetric organization of 
visual areas (that is, maps whose values reflect the orientation of the visuotopic 
gradient vector relative to the vector that points ‘geodesically’ towards V1, see 
Supplementary Methods 6.5). 

The classifier analyses were conducted using a standard machine learning train/ 
test/validation approach. The classifier was trained using the 210P subjects and 
tested against overfitting using the extra 29T subjects. The 210V subjects were 
used as the validation sample, and thus were not involved in the classifier training, 
testing, or the parcellation itself, and also shared no family relationships with the 
210P or the 29T groups. A short initial run of the classifier was used to identify 
features that the classifier was particularly sensitive to for each area (see below and 
Supplementary Methods 6.6). These features were compared in each individual 
subject with the group average pattern to exclude subjects that were potentially 
misaligned with the typical subject in this region (and hence for which the group 
defined training labels were likely inaccurate). This area-specific set of subjects 
in the 210P and 29T groups were excluded from the final classifier training of 
each area. The classifier’s output (ranging from 0 to 1) represents the likelihood 
that a given vertex in a subject is part of the area being classified or part of the 
surrounding cortex of the searchlight. Once the classifier training weights have 
been generated, it is possible to classify any subject who has the 112 multi-modal 
maps computed, including those whose areas are misaligned with the group 
(see Supplementary Results and Discussion 1.4). 

The trained classifier was applied to the 449 subjects and 27 repeat subjects 
to generate individual subject likelihood maps for each of the 180 areas in each 
hemisphere. These probability maps were combined by finding the largest 
probability for each vertex and then regularized within local neighbourhoods 
(see Supplementary Methods 6.7) to make an individual subject ‘winner-take- 
all parcellation. An area was considered to have been detected in a subject for 
the purposes of the areal detection measures (the overall classifier areal detection 
rate and the maps of areal detection rate for each area) if its size was between 
1/3 and 3 times the size of the original population-based parcel (a pragmatic 
threshold chosen prior to performing the analysis that tolerates modestly greater 
neuroanatomical variability across subjects than the empirical range reported in 
cytoarchitectonic studies*!”). Probabilistic maps of each area were then created 
by separately averaging the individual subject winner-take-all parcellation areas 
for the 210P and 210V subject groups. A group maximum probability map (MPM) 
parcellation was then created by assigning the identity of the maximum areal 
probability to each vertex. The reproducibility of the parcellation was assessed 
by correlating these two MPM maps and by computing a Dice coefficient. 
In both cases the parcellation was first turned into 180 concatenated binary 
ROIs per hemisphere (each area was represented by a separate map, ~30,000 
vertices per hemisphere, with ones for all vertices inside the area and zeros 
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for all vertices outside). The reproducibility of the individual subject hard 
parcellation maps was assessed similarly. For more details, see Supplementary 
Methods 7.2. 

Multi-modal areal fingerprints learned by the classifier were visualized using 
a classifier sensitivity metric. This metric was the partial derivative with respect 
to each feature of each area multiplied by the gradient magnitude of the feature 
(see Supplementary Methods 6.8). The measure indicates which areal features 
the classifier finds most informative when classifying a given area and whether 
increases or decreases in the value of the feature make the area more likely to be 
present. The sensitivity metric can be visualized both at the dense (vertex-wise) 
level for each feature and each area, or summarized at a parcel level. For each 
feature, the sensitivity metric was summarized at the parcel level by taking the 
maximum absolute value of the metric (finding the border where the feature was 
most influential) and using this maximum to represent the area in a parcellated 
or a matrix view, as shown in Supplementary Fig. 12 of Supplementary Results 
and Discussion 1.6. 
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Data reporting. No statistical methods were used to predetermine sample size. The 
experiments were not randomized. The investigators were not blinded to allocation 
during experiments and outcome assessment. 
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